function A409()
format long;
% Gaussian混合模型生成
% 假设所有数据（无论是否有标记）都由GMM生成，且每个类别对应一个Gaussian混合成分。

% ========== 图 1：二维 GMM 概率密度图 ========== %
mu = [1 2; -3 -5];                                                           % 定义了 二维高斯分布的均值向量
sigma = cat(3, [2 .5], [1 1]);                                               % 定义了 二维高斯分布的协方差矩阵，将两个 2x2 的协方差矩阵连接在一起，形成一个三维数组
% 第一个协方差矩阵 [2, 0.5; 0.5, 0.5] 对应第一个高斯分布。
% 第二个协方差矩阵 [1, 0; 0, 1] 对应第二个高斯分布。
gm = gmdistribution(mu, sigma);
subplot(2, 3, 1)
gmPDF = @(x, y) arrayfun(@(x0, y0) pdf(gm, [x0 y0]), x, y);
fsurf(gmPDF, [-10 10])
title('包含两个高斯分布的二维高斯混合模型（GMM）的概率密度')
xlabel('第一维度')                                                           % 这个维度表示数据在二维空间中的 水平坐标轴，即数据的第一维度。
ylabel('第二维度')                                                           % 这个维度表示数据在二维空间中的 垂直坐标轴，即数据的第二维度。

% ========== 图 2：GMM 拟合散点 + 等值线 ========== %
mu1 = [1 2]; Sigma1 = [2 0; 0 0.5];
mu2 = [-3 -5]; Sigma2 = [1 0;0 1];
rng('default');
X = [mvnrnd(mu1, Sigma1, 1000); mvnrnd(mu2, Sigma2, 1000)];
GMModel = fitgmdist(X, 2);
subplot(2, 3, 2)
y = [zeros(1000, 1); ones(1000, 1)];
h = gscatter(X(:, 1), X(:, 2), y);
hold on
gmPDF = @(x, y) arrayfun(@(x0, y0) pdf(GMModel, [x0 y0]), x, y);
g = gca;
fcontour(gmPDF, [g.XLim g.YLim])
title('GMM 在二维空间中对两个高斯分布数据的拟合结果')
xlabel('第一维度')                                                           % 这个维度表示数据在二维空间中的 水平坐标轴，即数据的第一维度。
ylabel('第二维度')                                                           % 这个维度表示数据在二维空间中的 垂直坐标轴，即数据的第二维度。
legend(h, '第一个高斯分布', '第二个高斯分布')
hold off

% ========== 生成 GPR 数据并 PCA ========== %
[gprData, gprLabels] = generateGPRdata();
classes = unique(gprLabels);
[coeff, score] = pca(gprData, 'NumComponents', 2);
disp('主成分系数矩阵（coeff）：')
disp(coeff)

% ========== 图 3：生成数据的原始分布 ========== %
subplot(2, 3, 3)
gscatter(gprData(:,1), gprData(:,2), gprLabels)
title('生成的 GPR 数据原始分布')
xlabel('反射振幅')
ylabel('介电常数')

% ========== 图 4：主成分 1 和主成分 2 的分布图 ==========  %
subplot(2, 3, 4)
gscatter(score(:, 1), score(:, 2), gprLabels)
title('PCA降维后（2个成分）的数据分布')
xlabel('第一主成分')
ylabel('第二主成分')

% ========== 图 5：GMM 拟合 PCA 得分 (k=2) ========== %
GMModel1 = fitgmdist(score, 2);
subplot(2, 3, 5)
h1 = gscatter(score(:, 1), score(:, 2), gprLabels);
h = gca;
hold on
gmPDF = @(x, y) arrayfun(@(x0, y0) pdf(GMModel1, [x0 y0]), x, y);
fcontour(gmPDF, [h.XLim h.YLim], 'MeshDensity', 100)
title('PCA降维后（2个成分）GMM拟合')
xlabel('第一主成分')
ylabel('第二主成分')
legend(h1)
hold off

% ========== 图 6：随机初始值拟合 ========== %
X2 = gprData(:, 1:2);
rng(10);
GMModel2 = fitgmdist(X2, 2);
subplot(2, 3, 6)
h = gscatter(X2(:, 1), X2(:, 2), gprLabels, [], 'o', 4);
haxis = gca;
xlim_val = haxis.XLim;
ylim_val = haxis.YLim;
d = (max([xlim_val ylim_val]) - min([xlim_val ylim_val])) / 1000;
[X1Grid, X2Grid] = meshgrid(xlim_val(1):d:xlim_val(2), ylim_val(1):d:ylim_val(2));
hold on
contour(X1Grid, X2Grid, reshape(pdf(GMModel2, [X1Grid(:) X2Grid(:)]), ...
    size(X1Grid, 1), size(X1Grid, 2)), 20)
uistack(h, 'top')
title('用随机参数初始化拟合生成 GPR 数据特征 1、2')
xlabel('特征 1 反射振幅')
ylabel('特征 2 介电常数')
hold off

% 显示估计的分量均值
muTable = table(GMModel2.mu(:,1), GMModel2.mu(:,2), 'VariableNames', {'主频', '衰减系数'});
disp(muTable)

% ========= 模拟的2分类GPR数据 ========= %
% 每个类别的样本数据由4个特征组成，特征由均值和标准差不同的正态分布生成
    function [gprData, gprLabels] = generateGPRdata()
        numSamplesPerClass = 50;
        numFeatures = 4;
        numClasses = 2;
        gprData = zeros(numSamplesPerClass * numClasses, numFeatures);
        gprLabels = cell(numSamplesPerClass * numClasses, 1);
        classNames = {'metal pipe', 'rock layer'};

        % 金属管
        gprData(1:numSamplesPerClass, :) = [... 
            normrnd(8.0, 1.0, [numSamplesPerClass, 1]), ...                  % 特征1反射振幅
            normrnd(1.1, 0.2, [numSamplesPerClass, 1]), ...                  % 特征2介电常数
            normrnd(120, 20, [numSamplesPerClass, 1]), ...                   % 特征3频率响应
            normrnd(0.1, 0.05, [numSamplesPerClass, 1])];                    % 特征4衰减系数
        gprLabels(1:numSamplesPerClass) = repmat(classNames(1), numSamplesPerClass, 1);

        % 岩层
        gprData(numSamplesPerClass+1:end, :) = [... 
            normrnd(3.0, 0.8, [numSamplesPerClass, 1]), ...
            normrnd(6.0, 1.0, [numSamplesPerClass, 1]), ...
            normrnd(40, 10, [numSamplesPerClass, 1]), ...
            normrnd(0.8, 0.15, [numSamplesPerClass, 1])];
        gprLabels(numSamplesPerClass+1:end) = repmat(classNames(2), numSamplesPerClass, 1);

        % 添加噪声
        noiseLevel = 0.05;
        gprData = gprData + noiseLevel * randn(size(gprData));
    end
end
